subroutine DQCoeff(xlist, maxorder, c)
! 作者：何光辉,
! 邮箱：flamehe@163.com
    !xlist: input, one-dimensional real list, with length n
    !maxorder : input, interger
    !c : output, differential quadrature coefficient array, with dimensions (0:maxorder, n, n)
    implicit none
    integer(4), parameter :: iwp=SELECTED_REAL_KIND(32)    
    integer(4) :: maxorder, i, j, m, n
    real(kind=iwp) :: xlist(:)
    real(kind=iwp) :: c(0:,:,:)
    
    n = ubound(xlist,1)
    c = 0.0_iwp
 
    do i = 1, n
        do j = 1, n
            if (i /= j) then
                c(1,i,j) = Mp(xlist,i)/((xlist(i) - xlist(j))*Mp(xlist,j))
            end if
        end do
    end do
    
    do i = 1, n
        c(1,i,i) = -sum(c(1,i,1:i-1)) - sum(c(1, i, i+1:n))
    end do
    
    do m = 2, maxorder
        c(m,:,:) = matmul( c(1,:,:), c(m-1,:,:) )
    end do
    
    do i=1,n; c(0,i,i) = 1.0_iwp; end do
        
    contains
    
    function Mp(xlist, i)
        implicit none
        integer(4), parameter :: iwp=SELECTED_REAL_KIND(32)
        real(kind=iwp) :: Mp, xlist(:)
        integer(4) :: n, i, j
        n = ubound(xlist,1); Mp = 1.0_iwp
        do j = 1, i-1
           Mp = Mp * ( xlist(i) - xlist(j) )
        end do
        do j = i+1, n
            Mp = Mp * ( xlist(i) - xlist(j) )
        end do
    end function Mp    
end subroutine DQCoeff